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ATS SIMULTANEOUS AND TURNAROUND RANGING EXPERIMENTS 

John S. Watson 
Barbara H. Putney 


ABSTRACT 

This report explains the data reduction and spacecraft position de- 
termination used in conjunction with two ATS experiments — Trilateration 
and Turnaround Ranging - and describes in detail a multilateration program 
that is used for part of the data reduction process. The process described 
is for the determination of the inertial position of the satellite, and for for- 
mating input for related programs. In the trilateration procedure, a geo- 
metric determination of satellite position is made from near simultaneous 
range measurements made by three different tracking stations. Turna- 
round ranging involves two stations; one, the master station, transmits the 
signal to the satellite and the satellite retransmits the signal to the slave 
station which turns the signal around to the satellite which in turn retrans- 
mits the signal to the master station. The results of the satellite position 
computations using the multilateration program are compared to results of 
other position determination programs used at Goddard. All programs give 
nearly the same results which indicates that because of its simplicity and 
computational speed the trilateration technique is useful in obtaining space- 
craft positions for near synchronous satellites. 
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ATS SIMULTANEOUS AND TURNAROUND RANGING EXPERIMENTS 
I. INTRODUCTION 

The purpose of this document is to explain the data reduction used in 'con- 
junction with two ATS spacecraft position determination experiments - Trilatera- 
tion and Turnaround Ranging - and to describe in detail one of the programs 
used for these experiments. In the trilateration procedure, a geometric deter- 
mination of satellite position is made from near simultaneous range measure- 
ments made by three different tracking stations (see Figure 1). Turnaround 
ranging involves two stations; one, the master station transmits the signal to 
the satellite and the satellite retransmits the signal to the slave station which 
turns the signal around to the satellite which in turn retransmits the signal to 
the master station. 

The Trilateration experiment used near simultaneous measurements of range 
to obtain a position vector of the satellite. In reality, three stations took measure- 
ments during the same ten minute period of time and by fitting this data with three 
polynomials as a function of raw data time one is able to create groups of simul- 
taneous measurements from the polynomials. See page 18 for details of the 
ranging pattern. These groups of simultaneous measurements were then fed 
into a Hughes trilateration program - Trilatron 1 (Reference 1), and out of each 
group of measurements one position vector of the satellite was obtained. The 
advantage of the trilateration technique is that it does not require a sophisticated 
force model that is required by conventional orbit determination programs. The 
trilateration method requires that the three tracking stations be located so as 
to give a well determinated triangle, with the subsatellite point preferably being 
located in the triangle. 

Turnaround Ranging uses a master station and some slave stations. The 
signal leaves the master station, travels to the satellite, returns to a slave 
station and is turned around to the satellite and then back to the master station, 
where it is recorded. By using this scheme, the slave station equipment need 
only be capable of turning the signal around to the satellite. Because the master 
station is tracking the satellite in its normal mode (two way delay time), the 
slave station range can be extracted from the four way delay time. 

In addition to describing the Trilateration and Turnaround Ranging experi- 
ments and their results, this report describes a more general multilateration 
program to determine the position of a satellite in synchronous orbit. In this 
program the position of the satellite is computed by knowing the ranges from 
three or more observation points, and the times of the observations which should 
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be as near the same time as possible. To obtain a solution, the program begins 
by making an initial estimate of the inertial coordinates of the satellite. For 
the first observation time, this initial estimate is obtained by using the trilatera- 
tion technique given in Reference 2. The initial estimate for each succeeding 
observation point is taken from the results of the previous one. From the inertial 
coordinates of the tracking station location and of the satellite, the program com- 
putes a range for each observation point and compares these to the observed 
ranges. By using a Taylor series expansion and a least squares process a 
converged set of the inertial coordinates of the satellite is then determined. 

If there are only three observation points, Cramer’s rule may be used to solve 
the set of three simultaneous equations for the position of the satellite. In 
addition, the multilateration program referred to in this report may be used to 
obtain a solution if there are more than three observation points. 
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H. DATA TYPES 


Two types of data were used. The two way delay data is the normal kind 
from the ATSR ranging system. It measures the time it takes for a signal to 
travel from the station to the satellite and return to the ground. This data 
type and its measurement is given in references 3 and 4. The other data type 
is four way delay time or Turnaround Ranging. This measurement signal leaves 
the master station, travels to the satellite, returns to the slave station, is 
turned around to the satellite and finally is returned to the master station. The 
master station must have tracked in the two way mode and in the same time 
period as well, in order to extract the slave station’s range at this time from 
the turnaround range. 


Trilateration Experiment 

In order to be able to trilaterate in the sense of the Hughes trilateration 
program, it is necessary to obtain a series of overlapping range values to three 
different stations. The Hughes Trilateration program then creates a position 
and velocity vector from this information. For detailed information about this 
program see reference 1. One of the functions of the software package described 
in this document is to create the required format from the polynomials repre- 
senting the raw data. See Reference 4 for raw data format information. 


Turnaround Ranging 

From the four way delay data a range to the satellite must be determined. 

The only station that needs to be fully equipped is the master station. The slave 
station need only to be able to turn the signal around. The output of this software 
package is a preprocessed series of ranges to be input to the trilateration pro- 
gram, and a DODS observation tape to be input to an orbit determination program. 


in. ALGORITHM AND FORMULAS 

This section will describe the formulas used and how they are derived to 
process this data. The details of the two way range data type are described in 
reference 3 and it is suggested that the reader read that document. The actual 
programmed equations will be marked with a starred number. It is hoped that 
the intermediate steps will make their derivation clear. The major problem in 
interpreting this data is that the time on the raw data message is not the proper 
observation time tag. Iterating so as to converge on the proper time tag and 
relating this time tag to the raw data time is the one complication of the proce- 
dure described in this algorithm. 
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A smoothing program creates a Chebyshev polynomial up to 12th degree 
of raw data and raw data time. See reference 3, Section 6.3. The spans of data 
have to be organized and read into this program in the groups that one wishes 
to have fitted with polynomials. There are sorting programs to do this. The poly- 
nomial coefficients and other details for creating the DODS Observation Tape 
are passed to the multilateration program where the ranges and proper time 
tags are then determined. 

A time T R , is chosen to be satellite time. The relationship between the raw 
data time and the ground received time is shown in the following figures. In order 
to find the proper measurement to correspond to the satellite time T R , one needs 
the raw data time T D that corresponds (the measurement polynomial is a function 
of raw data time). The connection between the times can be achieved by iterating 
to find the proper raw data time and therefore proper measurement to corre- 
spond to the satellite time,T R . The convergence is very rapid and few iterations 
are necessary. 

Let 

T d = raw data time 

8 2 = two way delay time raw measurement 

D m = total two way delay time of the master station 
S 4 = four way delay time raw measurement 
D 4 = total four way delay time 
D s = total two way delay time of the slave station 
N am = the ambiguity number for the master station 
N as = the ambiguity number for the slave station 
A A = the size of the ambiguity in time 
T r = the proper data time tag 
Tj = time the signal left the ground 
T f = time the signal returned to the ground 
c = speed of light 
S M = master station delay 
S s = slave station delay 
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One first chooses sets of data which overlap in time. For both trilateration 
and turnaround ranging it is a requirement to have the ranges from all the sta- 
tions involved at the same time. 

1* T r = first overlap time 


Since the polynomials are in terms of T D time an estimate of the T D corre- 
sponding to T r needs to be made. We will first consider the 2 way delay measure- 
ment. 


Since (T) 


T = 


T i + T f 


© T I = t d - 

© iy = t d + s 2 



T = T - 

R *D 


N AM AA S 2 



T r + 


n am aa 


and since S 2 is unknown at this point the first approximation to T D is 


2 * 


T d = T r + 


n am AA 


from equation 
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3* Then obtain S 2 at T D from the Chebyshev polynomial. 


Compute from equation 



4* 


T' = T - 
a r a d 


N»„AA S 2 


Compare 


5* If |T R /T' 
T d should be 


- 1 1 > e,T R and T R ' have not converged. According to equation (IT) 


T d = Tr + 


N AM AA S 2 


and we approximated it by 


2 * 


T - T 

a d a r 



AA 

__ 


so let 


6 * 



and try again from equation 3*. 

5* If | T r / T r - l| < e T r and T R have converged, then compute 
7* D„ = S„ + N.„ A A - S,, the master station total delay. 

M 2 AM M 

If all the data is two way data one would repeat this process with each station. 

At this point, one would have the proper measurement from each station at 
the same time. One could solve for a position vector at time T R (trilaterate) if 
one has ranges from three stations. 

If some of the data is four way (turnaround) data the following algorithm 
must be used. 
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The four way data contains the sum of ranges to two stations. It is necessary 
to extract the range to the slave station from this data. In order to do this, one 
must have tracked the satellite from the master station alone at approximately 
the same time. This is only necessary because the orbit determination programs 
currently do not accept the summed range as a data type. Due to this restric- 
tion, one relies on the smoothing polynomial to fill in the gaps in the data. It is 
assumed at this point that the master station total delay is known at time T R , 
as obtained in the previous algorithm. It is now necessary to find the four way 
total delay, D 4 and then to compute the total two way delay for the slave station 

Choose the same T R as before. 

Since (IT) 

© 

© 

© 

© 

and since S 4 is unknown at this point the tirst approximation to t d is 

9* Then obtain § 4 at T D from the Chebyshev polynomial. 

Compare S 4 and S 2 . 

If S 4 1 b 2 then an extra ambiguity had to occur because of the combined 
effect of two stations. This means equation (lo) needs an extra AA/2 added to 
it and 


=T d -N am AA-N as AA 


T f = T d + § 4 


= T °-( 


^AM^A+N^AAX ^ 


T D “ Tr + 


N am A A+N as AA\ S 
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^ n «i.aa+n m &a + aa \ 


for this case. 
Therefore 


10 * 



+ 


AA 

2 


and since the initial T D had to be far off another 10* S 4 is obtained from the 
Chebyshev polynomial at the new T D . 

D 4 is then computed. 


11 * 


D 4 = 8 4+N am AA + N as AA + AA-S m -S s 


/N A A + N a<? A A + A A\ 8 

12* T r = T d _ ( 2 ) + ~T fr ° m ec l uation (11. 


If the new S 4 > § 2 recompute T R and D 4 as if that were the original situation 
as below. When S 4 > S 2 compute 


13* 




AA + n as 
2 




from equation 



and 14* 


D 4 = S4+(N am AA + N as AA)-S m -S s 


If 15* I T r /Tr - 1 1 1 £, T r and T R have converged then 
16* Ds = D 4 - D m 


If 15* | T r /T^ - 1 | > e , T r and 


have not converged then 


let 16* 


= T - 


2 


and try again from equation 3*. 
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After two way delays are determined for each station the range R for each 
station is computed 

/ 

17* R = ~2 D. where D = D M or D g 

The output is properly formatted and the data is time tagged, thus completing 
this task. 

IV. COMPUTATION OF THE INERTIAL COORDINATES OF 
THE TRACKING STATIONS 

From the results of the previous section, we now have a set of ranges (three 
or more) at the same time. The function of the multilateration program is to 
compute the sub-satellite point at this time. 

> j 

The multilateration method takes the set of three or more simultaneous 
ranges, uses a least squares solution to compute the inertial satellite position, 
and then computes the sub-satellite point in the earth fixed system. These 
items are computed according to the formulas presented in sections IV, V, and 
VI. 

However, first one needs to know the coordinates of the ranging stations. 

The method for computing the inertial coordinates of the tracking stations is 
the one presented on pages 16-18 of Reference 5. 

The following data are needed to compute the coordinates of the tracking 
stations: 

(\ 0 ) = the hour angle of the first point of Aries 

(\ E ) = the geodetic longitude of the terrestrial tracking station in radians, 
as measured eastward from Greenwich (a negative sign must be 
prefixed if measured westward from Greenwich) 

(# D ) = the geodetic latitude of the station in radians, measured as positive 
north of the Equator, and as negative south of the Equator 

(H' ) = the altitude of the station in feet, measured positive above sea level 
and negative below sea level 

(co) = the angular velocity of rotation of the earth in radians per hour 

(A T) = the difference in hours between the observation time and midnight 
preceding the observation time 
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(f) = the flattening coefficient of the earth 
(e c ) = the geocentric latitude of the tracking station in radians 

(P) ~ the geocentric distance of the tracking station, in units of earth 
radii 

(e) = the eccentricity of the earth 

The geodetic longitude (A- E ), geodetic latitude ($ D ) , and height (H) of the 
tracking station are known. From these the inertial geocentric coordinates of 
the tracking station in spherical coordinates Q, 0 , s ) can be computed. These 
spherical coordinates can then be converted to a Cartesian system of coordi- 
nates (x T , y T , z T ). 

In the meridian section of the earth through an observer, the position of the 
latter relative to the earth’s center can be expressed in rectangular coordinates 
as: 

p sini9 G = Ssin0 D , p cos 0 Q = C cos 0 D . 

These serve to define the auxiliary functions S and C. 

C = [cos 2 6 > d + (1 _ f) 2 sin 2 0 D ] -1/2 

S = ( 1 - f ) 2 C 

H = (4.77865 x 10~ 8 ) H' (converts feet to earth radii) 

s ° = arctan [(fri)] tan 8 ° 

P = [(S + H) 2 sin 2 0 D + (C + H) 2 cos 2 0 D ] 1/2 
5 =\ 0 +U (AT) +\ E 

(S) = the angle in radians between the vernal equinox and the 
observation meridian plane 
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X T = p cos d Q cos S 
yj = £ cos # G s in 8 


z T = p sin d Q 

(x T , y T , z T ) = the geocentric coordinates of the tracking station in units 
of earth radii 

We now have the station locations in quantities and units that we can use in 
the program. 


V. COMPUTATION OF THE INERTIAL COORDINATES 
OF THE SATELLITE 

This multilateration method computes the inertial coordinates of the satel- 
lite by a least squares iteration procedure, given the ranges from the satellite 
to three or more tracking stations, the inertial coordinates of the tracking 
stations, and an initial estimate of the inertial coordinates of the satellite. This 
section describes the mathematical method for computing the inertial coor dinat es 
of the satellite. We refer to this method as the multilateration method. 

(R 0 ) = observed range in earth radii 
(R c ) = computed range in earth radii 

(x, y, z) = the present estimate of the inertial coordinates of the satellite, 

obtained initially by using the trilateration method given in Refer- 
ence 2. 

(x’, y’, z') = the new estimate of (x, y, z) 

Rc = ( x - X T ) 2 + (y - y T ) 2 + ( z - Z T ) 2 

A function of three variables may be expanded in a series by Taylor’s formula 
in the form: 


f (x + Ax , y+Ay, z + Az) - f(x, y, z) 


A 3f A 3f A 

= — Ax + — Ay + — Az 
ox oy dz 


+ higher order terms 
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or Af=^-Ax+^Ay + ^Az, ignoring the higher order terms 

ox oy oz 


since 


R c - f ( x > y. z ) 


we can write 


, BR 2 BR 2 

A R 2 = — — A x + —— A y 
ox oy 


BR 2 

+ — -A Z 
oz 


where 


AR 2 = R 2 q - R 2 c 


3R 2 

— — = 2x - 2x 
ox 1 


BR 2 

V - = 2y " 2Yt 


BR 2 c 


= 2z - 2z, 


Therefore 


or 


AR 2 = (2x - 2x T ) Ax 

+ (2y - 2y x ) Ay + (2z - 2z T ) Az 
\ AR 2 = (x - X T ) Ax + (y - y T ) Ay 


+ (z - z_) Az 


(5.1) 


A least squares routine is then used to compute Ax , Ay , Az. The least 
squares method used is the one described on pages 61-66 of Reference 5. 
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Then 


x' = x + Ax 
y' = y + Ay 
z' - z + Az 


The program then returns to Equation 5.1 with x, y, z being replaced by x‘, 
y’, and z'. AR 2 is re-evaluated. This iterative procedure is continued until the 
standard deviation of fit of 1/2 AR 2 reaches the desired tolerance, and we thus 
have the inertial position of the spacecraft. 


VI. COMPUTATION OF THE LONGITUDE AND LATITUDE 
OF THE SUB-SATELLITE POINT 

The method for transforming the inertial position coordinates (x, y, z) of 
the satellite, as computed above, to the East longitude (k s ), the geodetic latitude 
(<9 d ), and the elevation of the satellite (H s ) above and normal to the adopted 
ellipsoid, and the geocentric radius magnitude ( p ) is given in this section. 


The symbols used in this section are the same as those used in Section IV, 
except they refer to the sub-satellite point instead of the tracking station location. 
The method used in this section may be found in Reference 1. 


The converged values of the geocentric coordinates of the satellite (x, y, z) 
have been found using Section V. Proceeding from these knowns, we compute 
the unknown longitude, latitude, and elevation as follows: 


y0 2 = x 2 +y 2 +z 2 


e 2 = 2f - f 2 
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cos & Q 


(X 2 + y 2 ) 

_ P 


sin & G - 


z 

P 


tan 0 Q 


sin <9 q 
cos 9 q 


tan 9 q = (1 - f) 2 tan 0 D 


d Q = arc tan 


tan 9 Q 
(1 - f) 2 


cos S = 


p cos 9 C 


s in 


S = 


p cos 9. 


X s = 8 - \ Q -"(AT) 


if X. is greater than 180°, the program will take X s as 


\ s = \ s - 360‘ 


A0 g = O 


0 o = e o-^ e o f 6 ' 1 ) 


1 - e 


1 1/2 


Ll - e 2 cos 2 9 G -i 
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t 

If A(9 g has not stopped varying return to Equation (6.1). 

H s ~ 3 e H 

(k s ) = the geodetic longitude of the sub-satellite point (+ indicates east 
of Greenwich and - west of Greenwich) 

(<9 d ) = the geodetic latitude of the sub-satellite point (+ indicates north 
of the Equator and - south of the Equator) 

( a e ) = radius of the earth in kilometers 

(H s ) = the height of the satellite above and normal to the adopted 
ellipsoid. 

Thus we have computed the longitude and latitude of the sub -satellite point, 
and the height of the satellite. 

VII. RESULTS 

Simulated Data 

The multilateration program described herein has been tested using simu- 
lated data for the ATS-1 satellite. The epoch used was the 4th of April 1969 at 
0 hours, 0 minutes, and 0 seconds. 

Some of the orbital parameters were: 

semi-major axis = 42,166.5 km. 
eccentricity = .000229 
inclination = 1.4° 

period = 23.94 hr. 
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The test cases listed below were run using the three observations given 
below. 


Observation Time 


Tracking Station 

yr. 

mo. day 

hr. 

min. 

sec. 

Range 

Mojave, Calif. 

69 

4 

4 

0 

37 

0 

38,133.6887 km, 

Rosman, N. C. 

69 

4 

4 

0 

37 

0 

40,633.4199 

Toowoomba, Australia 

69 

4 

4 

0 

37 

0 

39,518.6607 


The true values of x, y, z as given by the program which simulated the data 
were: 

x = 26,297.9477 km. 
y = 32,944.5272 
z = -573.1739 


The initial values used for the three test cases were: 

x y z 

Case 1 +28,701.7425 km. +28,701.7425 km. -4464.7155 km. 

Case 2 +19,134.4950 +22,323.5775 -6314.3833 

Case 3 +25,512.6600 +31,890.8250 -574.0348 

All the cases converged to the same values: 

x = 26,298.1470 km. 
y = 32,944.3690 
z = -573.1558 

If we compare these values to the given values, we get: 

Ax = -200 m. 

Ay = +158 
Az = -18 
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Where the difference is taken as the given value minus the computed value. 

The sub-satellite point is then found to be: 

Longitude = -141.66° 

Latitude = -0.784° 

Height of satellite = 35,779.31 km. 


The program gave similar results when Cramer’s method was used instead 
of a least squares procedure. For example, using Cramer’s method for Case 1 
gave the following results: 


X = 

26,298.1546 km, 

y = 

32,944.3582 

z 

-573.1764 

Ax = 

-207 m. 

Ay = 

+169 

Az -- 

+2 

Longitude = 

-141.660° 

Latitude = 

-0.784° 

Height of satellite = 

35,779.92 km. 


The largest error in the test cases was 207 meters. These test cases do 
not agree exactly because the program generating the simulated data and the 
multilateration program did not use the same sidereal time. These test cases 
were run before the trilateration method for computing the initial estimates of 
x, y, z had been incorporated into the multilateration program. Guided by these 
results, the program was then run for three cases using real satellite data. 

Real Data 


Case 4 

The observations for this case were taken on October 1, 1969 by three ATS 
ground stations: Rosman, North Carolina; Mojave, California; and Toowoomba, 
Australia. Mojave was the prime or two-way delay station; whereas, Toowoomba 
and Rosman were four-way delay stations. Two-way delay means that the signal 
is sent from the station to the satellite and is then returned to the station. 
Four-way delay means that the signal is sent from station A to the spacecraft, 
from the spacecraft to station B, from station B to the satellite, and from the 
satellite to station A. Station A would be the master station in this case. 
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The three stations tracked the spacecraft during a ten minute span of time 
as illustrated in Figure 2. The data rate was one observation per second and 
each station ranged for about one minute each time it tracked. 


© © © © © © © 



Rosmon 


© 


Toowoombo 


Figure 2. 


Some of the orbital parameters used are: 

semi-major axis = 42,168.1 km. 
eccentricity = .00229 
period = 23.94 hr. 

The initial values used for the inertial coordinates are: 

x = +5740.3485 km. 
y = -41,458.0725 
z = -3507.9907 

The results are then compared to the results obtained by the Refined World Map 
Program (Reference 6). 


Observation Time Satellite Position 


Program 

hr. 

min. 

sec. 

Long. (+E) 

Lat. (+N) 

HT 

ATS 

3 

47 

56.00 

-148.710° 

-.484° 

35782.0 

WMAP 

3 

48 

0.00 

-148.650 

-.489 

35787.0 

ATS 

3 

53 

56.01 

-148.734 

-.536 

35780.7 

WMAP 

3 

54 

0.00 

-148.655 

-.537 

35788.7 

WMAP 

4 

41 

0.00 

-148.660 

-.896 

35789.0 

ATS 

4 

41 

35.00 

-148.779 

-.822 

35787.0 

WMAP 

4 

42 

0.00 

-148.660 

-.903 

35790.0 


Here, ATS refers to the multilateration program developed in this report. 
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The longitude agrees to within .1°, the latitude within .01° and the height 
within 10 kilometers. Some of the error is probably due to the delays in the 
ground equipment at the master station not being measured accurate enough. 


Case 5 


The observations for this case were taken on December 23, 1969 by the ATS 
ground stations: Rosman, Mojave, and Toowoomba. This time the stations made 
simultaneous observations. All the stations were therefore two-way delay stations. 
They give an initial estimate for x, y, z of 

x = +21,934.0285 km. 
y = -35,990.6776 
z = -1103.2987 

Below is the comparison of the results of the multilateration program (called 
ATS in the charts) with the ORB-1 and WMAP programs (Reference 6). 


Observation Time 


Program yr. 

mo. 

day hr. 

min. 

sec. x 

x 

z 

ATS 

69 

12 

23 

0 

22 

0.02 25,677.7275 

km.-33, 421. 9067 km. 

-1124.1973 km, 

ORB-1 

69 

12 

23 

0 

22 

0.00 25,683.7708 

-33416.8086 

-1125.5546 

ATS 

69 

12 

23 

0 

23 

1.02 25,826,0626 

-33,307.2853 

-1128.8897 

ORB-1 

69 

12 

23 

0 

23 

0.00 25,829.6496 

-33,304.0089 

-1130.0527 


Observation Time Satellite Position 


Program 

yr. 

mo. 

day 

hr. 

min. 

sec. 

Long. (+E) 

Lat.(+N) 

Ht 

ATS 

69 

12 

23 

0 

22 

0.02 

-149.340° 

-1.538° 

35783.8 km, 

WMAP 

69 

12 

23 

0 

22 

0.00 

-149.327 

-1.531 

35784.0 

ATS 

69 

12 

23 

0 

23 

1.02 

-149.340 

-1.545 

35783.8 

WMAP 

69 

12 

23 

0 

23 

0.00 

-149.327 

-1.537 

35784.0 


The longitude agrees to within .1° the latitude within .01°, and the height 
within 1 kilometer. 
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Case 6 


The observations for this case were taken on February 20 and 21, 1970. 
The observations were taken simultaneously by the three tracking stations of 
Mojave, Kashima, and Toowomba. 


Below is a comparison of the multilate ration program (called ATS in the 
charts), with the TV-1, ORB-1, and WMAP programs. 


Observation Time Satellite Position 


Program 

yr. 

mo. 

day 

hr, min. sec. 

X 

x 

z 

ATS 

70 

2 

20 

4 

40 30.26 +148 43. 66km. 

+39476 .84km. -349.72km, 

TV-1 

70 

2 

20 

4 

40 30.26 +14844.79 

+39476.42 

-349.72 

ORB-1 

70 

2 

20 

4 

40 30.00 +14844.87 

+39476.33 

-350.27 

ATS 

70 

2 

20 

10 

40 29.26 -39487.29 

+14695.89 

+1642.31 

TV-1 

70 

2 

20 

10 

40 29.26 -39486.87 

+14702.83 

+1642.26 

ORB-1 

70 

2 

20 

10 

40 30.00 -39487.87 

+14699.94 

+1641.75 

ATS 

70 

2 

20 

17 

10 30.25 -9206.24 

-41136.44 

+120.00 

TV-1 

70 

2 

20 

17 

10 30.25 -9207.42 

-41136.18 

+120.00 

ORB-1 

70 

2 

20 

17 

10 30.00 -9207.20 

-41136.14 

+119.27 

ATS 

70 

2 

20 

21 

10 28.28 +31088.54 

-28439.97 

-1391.48 

TV-1 

70 

2 

20 

21 

10 29.29 +31087.26 

-28438.56 

-1391.55 

ORB-1 

70 

2 

20 

21 

10 30.00 +31096.58 

-28436.10 

-1392.39 





Observation Time 


Satellite Position 

Program 

yr- 

mo. 

day 

hr. 

min. sec. 

Long. (+E) 

Lat. (+N) 

Ht. 

ATS 

70 

2 

20 

4 

40 30.26 

-150.44° 

O 

00 

• 

1 

35798.6 km. 

TV-1 

70 

2 

20 

4 

40 30.26 

-150.44 

- .48 

35798.6 

WMAP 

70 

2 

20 

4 

40 30.00 

-150.43 

- .47 

35798.4 

ATS 

70 

2 

20 

10 

40 29.26 

-150.49 

+2.23 

35789.2 

TV-1 

70 

2 

20 

10 

40 29.26 

-150.49 

+2.23 

35789.2 

WMAP 

70 

2 

20 

10 

40 30.00 

-150.49 

+2.23 

35789.0 

ATS 

70 

2 

20 

17 

10 30,25 

-150.46 

+0.16 

35776.0 

TV-1 

70 

2 

20 

17 

10 30.25 

-150.46 

+0.16 

35776.0 

WMAP 

70 

2 

20 

17 

10 30.00 

-150.45 

+0.16 

35775.8 

ATS 

70 

2 

20 

21 

10 28.28 

-150.45 

-1.89 

35778.5 

TV-1 

70 

2 

20 

21 

10 29.29 

-150.45 

-1.89 

35778.5 

WMAP 

70 

2 

20 

21 

10 30.00 

-150.45 

-1.89 

35778.4 
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The longitude agrees to within 0.1°, the latitude within .01°, and the height 
within 1 kilometer. The TV-1 and the multilateration programs did not give the 
same results, because the TV-1 programs used mean sidereal time and the 
multilateration program used apparent sidereal time. 

From these results it can be seen that there is excellent agreement between 
the various position determination methods. The trilateration method has ad- 
vantages over conventional orbit determination techniques in that it is a geo- 
metric solution and does not require a sophisticated force model, nor does it 
require a number of iterations to obtain an orbit as do conventional orbit de- 
termination techniques. Each position vector takes less than five seconds CPU 
time on the 360/91 when using the trilateration program. The time required by 
conventional orbit determination programs depends on the length of the arc 
needed to determine the orbit, the number of iterations necessary to converge, 
and the force model used. Some disadvantages of the described trilateration 
technique are; an orbit is not determined and therefore orbit prediction can not 
be made, and there must be mutual visibility by three stations capable of track- 
ing almost simultaneously the particular satellite. For additional information or 
the future development of the trilateration technique the reader should refer to 
references 8 and 9. 


VHI. CONCLUSION 

From the results of the cases given in the previous section, it can be seen 
that the multilateration programqgives good results. The Multilateration, TV-1, 
ORB-1 and WMAP programs compared very favorably. The longitude always 
agreed to within .1 degrees, and the latitude always agreed to within .01 degrees. 
The height agreed to within 1 kilometer, except on the turnaround ranging case. 

A possible explanation for this discrepancy is that the measurements of the time 
delay in the ground equipment at the slave station is inaccurate. The results 
further indicate that spacecraft positions for synchronous satellites can be ob- 
tained using less computer time and less computer memory if one uses the 
multilateratinn instead of using conventional orbit determination techniques. 
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APPENDIX 
Program Information 


A. Operating Instructions for Multilate ration Program 

Below is an explanation of the data cards, with some typical values in 
parenthesis 

1. Card #1-5 variables in the format (5E15.8) 

a. f = the flattening coefficient of the earth (.33528919) 

b. TOL = the tolerance that is used to determine when the least squares 

program has converged to a solution (.1 x 10“ 8 ) 

c. DELA = the size of the ambiguity in time in seconds (.125) 

d. REARTH = the radius of the earth in kilometers (6378.165) 

e. C = the speed of light in km./sec. (2.997925 x 10 s ) 

2. Card #2-7 variables in the format (2(15, 5X), F5.1, 15, 5X, E15.8, 215) 

a. NOOBS = the number of tracking stations being used in this run. The 

program with slight modification will handle up to 10 stations. 

(3) 

b. NOPRST = the station number of the prime or two way delay station. 

If all stations are two way delay stations this number 
would be set to -1. (47) 

c. DELTAT = the time increment in minutes at which the data points will 

be printed out in each time interval. (1) 

d. NOPTS = the number of data points that are to be computed in each 

time interval. This option would override the value of 
DELTAT. This option is usually not used and the value is 
set to -1. (-1) 
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e. EPSLON = the tolerance which is used to determine when the program 

has converged on the proper time (1.0 x 10 ~ 18 ). 

f. IMON = the month in which the observations begin (10). 

g. IDAY = the day of the month on which the observations begin (1) 

3. Card #3 to card # (2 + NOOBS) - 3 variables in the format (215, 5X, 

E15.8)) 

a. NUMST = the number of the tracking station (47). 

b. NA = the number of ambiguities (2) 

c. STNDEL = the station delay time in micro-seconds (1.75). 

4. Card #(2 + NOOBS) to Card # (32 + NOOBS) - 1 variable in the format 
(7X, E15.8). 30 cards total 

a. ALAMDM = the hour angle of the first point of Aries at midnight 

Greenwich mean time in radians. The first card in this 
set contains the hour angle for the IMON and IDAY on 
card #2. Each subsequent card has the hour angle for the 
next day. As the program is presently set up 30 cards are 
required, however if there is data only on one day only 
the first card need be filled in the other 20 may be blank. 

If NOOBS on card #2 is 2 or less all 30 cards may be 
blank. (.16668687, value for Oct. 1, 1969, obtained 
from Reference 8.) 

Last 4 cards give the coordinates of the four ATS tracking stations in the 
format (13X, E15.8, IX, E15.8, IX, F8.1). 

a. ALAMDE = the geodetic longitude of the tracking station in radians 

as measured positive eastward from Greenwich 
(+4.2430894) 

b. THETAD = the geodetic latitude of the tracking station in radians as 

measured positive north of the equator (+.61665814). 

c. ALT = the altitude of the tracking station in feet measured as positive 

above sea level (+3072.8). 

Note: The name of the station is given in columns 3-8 (Mojave in the 

example given in a thru c above), but is not used in the program. 
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B. Program Flow Chart 

The program is presently set up for the ATS satellites and must 
have the four tracking stations in the order given below. 

1. Kashima, Japan 

2. Mojave, California 

3. Rosman, North Carolina 

4. Toowoomba, Australia 

If any other stations but these four stations are used the multilateration 
program has to be slightly modified, because the station numbers are hard 
coded in the program. 
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C. Multilateration Program 

Thplicit REAL*8(A'-H,0-Z,*) 

REAL MAXR,MAXRO*MAXXY»MAX,MIN 

INTEGER CC1,CC2,CC3,CC4,C1SAV,C2SAV,C3SAV,CASAV,U,V 
RE AL*8 OTIME,STIO,OeSV,TRANOM,FICRP,OBCV 
INTEGER *4 ISATN.IOBN 

REAL *4 08WT,C0C,P0C,SDFT,FNRRES, TMER, TCR V.ROBtoT , RROBWT , A08WT 
INTEGER*2 I ESN, IDRPN, IOT,IEF,ITTI,ITCI,IOCI 

DIMENSION A ( 10,13), ACI13), ALAMCNI30), CLMQH10.3), COEFF(I3), 
lOELAY(LO), DELR (10, 1 ) , CLTAQH3.4), END(IO), NDAYI12), NDEG(IO), 
2NCMST(10), R (90 ) , RCBS(IO), STACCN(4,3), START(IC), TYPE(IO), 
3WTI1C), XT ( 1C) , YT(10), ZT(IO), RTVH200, 3), TVTIMEI200, 3), 

A 1 1 1 ( 2 ) , STATNM(IO), STOFIT(IO), toTOB(lO), DELOUT ( 2C0 , 10 ) , 
5NSTATI10), STOR (13), ST ADEL < 100 ) , DELD0D1200, 10) , NAMB(IOO) 
DIMENSION NAMBPRI 10) 

EQUIVALENCE (FIDRP, 111) 

NAMELIST/NAMl/N, NUMST, Y, NAM8, DELA, STNDEL. DEL 2 to , D2toAY, DELAY 
1 , Y 2 , Y4 

NAMELIST/NAM2/N, NUMST, Y, NAM8 , DELA, STNDEL, DEL2W, D2w AY , DELAY 
1 , Y2, YA 

NAMELI ST/NAM3/N, NUMST, Y, NAMB, DELA, STNDEL, DEL 2W , D2WAY, DELAY 
1 , Y2, YA 

C CORRECTION FOR LIGHT TIME AND AMBIGUITY REQUESTED 

INTEGERS IAMBRS 
DIMENSION I AMBR S ( 3 ) 

DATA IAMBRS/ 16, 22,3 1/ 

C RANGE RR ANGLES CORRECTICNS TO BE MADE 

DATA TRANOM, FIDRP ,OBS V , CB toT ,C CC , PCD , F NRRE S , TME R , TC R V , I E SN , I DR PN , 
1I0T, IEF, ITT 1, I TCI , I OC I/S*C. ,7*0/ 

196 IOT = 1 

AE = 6378.165 

R ADC EG = 57.2957795 

NDAY(l) = 31 

NDAYI2) = 28 

NO AY ( 3) =31 

NDAY(A) = 30 

NOA Y ( 5 ) = 31 

NDA Y ( 6 ) = 30 

ND A Y ( 7 ) = 31 

NO A Y ( 8 ) = 31 

NOA Y ( 9 ) = 30 

NOA Y ( 10 ) = 31 

NDA Y (11) = 30 

NOA Y ( 12 ) = 31 


IF NCPTS = 1 OR > l, .THE PROGRAM toILL COMPUTE A RANGE FOR THE 
NOPTS SPECIFIED. NO CCMfUTATI CN CF INERTIAL COORDINATES OR 
OF SUB-SATELLITE POINT WILL BE MAOE IN THIS CASE. 

NCPRST = -1 FOR THIS CASE. SET NOPTS = 0 FOR REGULAR CASE. 

IN THIS CASE THE PROGRAM toILL COMPUTE THE TIME CCMMCN TG ALL 
STATIONS AND COMPUTE RANGE, AND SUB-SATELLITE POINT FOR THIS TIME 
AND CONTINUE TO COMPUTE AT SPECIFIED INTEGRALS OF TIME. 

REAC ( 5 , 199 ) F, TOL, DELA, REARTH, C, NCObS, NOPRST, DELTA!, NOPTS, 
1 EPSLCN, I MCN , I DAY 

199 FORMAT! 5E15.8/2I 15, 5X) , F5.1, 15, 5X, E15.8, 215) 

WR I T E ( 6 , 5000 ) F, TOL, DELA, N008S, NOPRST, DELMIN, REARTH, C 
WRITE(6,5019) 

5C19 FORMAT ( IHO, 37X, 17HTRACK ING STATIONS) 

DO 195 I = 1, NOOBS 
RE AC ( 5 , 194) NUMST, NA, STNCEL 
19 A F ORM AT (215, 5X, E15.8) 

STNCEL = STNDEL ♦ 1. CD-06 
WRITE(6,5020) NUMST, NA, STNDEL 
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15, 5X 


5C20 FORMAT ( 1H0 , _ 

117H NO. OF STATION =, 15, 5X, 1 AH NO. OF AMB. = 

216H STATION CELAY =, E2C.8, 2X, AHSEC . // ) 

STACEL(NUMST) = STNCEL 
195 NAMB(NUMST) = NA 
DELM IN = DELTAT 

IF NOPRST = -1, ALL STATIONS ARE TWO WAY DELAY STATIONS 

RE AO ( 5 « 200 ) (ALAMOM(t), 1=1, 30) 

2C0 FORMAT ( 7X, E15. 8) 

REAC IN TRACKING STATION LOCATIONS 

201 FORMAT ( 13X , E15.8, IX, E15.8, IX, F8.1) 

REAC (5,201) <(STACCN(I,J), J = 1,3), 1 = l.A> 

I OBN = 0 

TEMP2 = OELA / 86A0C.C 

DELTAT = OELTAT / 1AA0 .0 

N = 1 

NOTIM = 0 

NOITER = 1 

III = 0 

NOREAD = 0 

JJ =1 

IF(NOPTS) 3, 3, A 
3 NK = NOOBS 

GO TO 5 
A NK =1 

5 00 10 1 = 1, NK 
NOREAD = NCREAD *1 
INDIC = -1 

REAC l 19 , END=8998 ) STTIM, FINTIM, NUMST, OATYP, NDGREE, COEFF, 

1 1 SAT , ROBWT, SOFT, STNANE , FIORP, IOCI 

THIS MEANS NO GATING IS DONE BY OODS, NOTICE THIS IS GOOD FOR 
RANGE ONLY, 00 NOT USE IF RANGE RATE IS ALSO USED 

111 ( 2 ) = 11112 ) ♦ 2 

NOMST ( I ) = NUMST 

IF ( NCREAD - NOOBS) t, 6, 7 

6 III = III ♦ l 
NST AT (III) = NOMSTI I) 

7 START! I ) = STTIM 

END ( I ) = FINTIM 

TYPE ( I ) = CATYP 

NDEG(I) = NDGREE 
STATNM ( I) = STNAME 
STDFIT(I) = SOFT 
HTOBU) = ROBWT 

KK = NDEGl I) + 1 

DO 11 J = 1, KK 

11 A( I ,J) = COEFF ( J ) 

12 I F ( J - 13) 13, 10, 10 

13 J = J + 1 

A ( I , J S =0. 

10 WR IT E ( 6 « 500 1 ) I, START(I), ENO(I), NOMST(I), TYPE(I), NCEG(I), 
1 { A 1 1 , J) , J = l, 13) 
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5C01 FORMAT! 1HO, 31 
110H START 
210H TYPE 
31 OH A1 2 ) 

410H A( 5 ) = 

510H A<8) 

610H A(U) 

GO TO 215 
139 INOCT = +1 


I = 13// 
E20.8//10F END 
E20.8//10F NOEG 
E20.8//10F A( 31 
E20.8//10F A( 6) 
E20.8//10F AT 9 ) 
E20.8//10F A ( 12 1 


= E20.8//10H NOMST 
= 15 //10H All) 

= E20.8//10H A C A ) 

= E20.8//10H At?) 

= E20.8//10H At 10) 
* E20.8//10H A ( 13 ) 


IF(NGPTS) 415, 415, 140 

140 M =1 

ENDTIM = FINTIM 

IFtNCPTS - 2) 141, 142, 143 . 

14 1 T = STTIM ♦ < <ENDTIF - STTIM) / 2.) 

DELTAT « ENDTIM 

GO TG 25 

142 DELTAT = (FINTIM - STTIM) / 3. 

DELTAT = DELTAT - TEMP 3 

T = STTIM ♦ DELTAT 

ENDTIM = FINTIM - (OELTAT / 2.) 

GO TO 25 

143 DELTAT = (ENDTIM - STTIM) / (NOPTS - 1.) 

OELTAT = DELTAT - TEMP3 

T = STTIM 

GO TO 25 

415 BEGTIM = START ( 1 ) 

ENDT IM = END! 1 ) 

DO 19 I = 2, NOOBS 
IF<START(I) - BEGTIM) 17, 17, 16 

16 BEGTIM = START ( I ) 

17 IF* ENC (II - ENOTIM) 18, 19, 19 

18 ENDTIM = ENO( I ) 

19 CONTINUE 

IFtBEGTIM - ENOTIM) 490, 490, 405 


= 15// 

= E20.8// 

= E20.8// 

= E20.8// 

= E20.8// 

= E20. 8/ / ) 


* 


NO OVERLAP PERIOD FOUND 


405 MM = NCOBS ♦ 1 

INDIC = ♦ l 

NOREAD = NOREAO ♦ l 

REAC(19,EN0 = 8998) STTIM, FINTIM, NUMST, DATYP , NDGREE , COEFF, 
IISAT, ROBWT , SOFT, STNAME , FICRP, I OC I 

THIS MEANS NO GATING IS DONE BY CODS, NOTICE THIS IS GOOD FOR 
RANGE ONLY, DO NOT USE IF RANGE RATE IS ALSO USED 

111 ( 2 ) = 111 ( 2 ) ♦ 2 

START(MM) = STTIM 
END ( MM ) = FINTIM 

TYPE(MM) = CAT YP 
NOMST ( MM ) = NUMST 
STATNMl MM) = STNAME 
NCEG(MM) = NDGREE 
STOFIT(MM)= SOFT 
WTOB(MM) = ROBWT 
KK = NDEG(MM) ♦ 1 

DO 451 J = 1, KK 

451 A ( MM , J ) = COEFF(J) 

452 I F ( J - 13) 453, 450, 45C 

453 J = J ♦ 1 

A ( MM , J ) = 0. 

450 WRITE! 6 ,500 1 ) MM, S TART (MM ) , END ( MM ) , NOMS TIMM), TYPEIMM), 
INDEG(MM), ( A ( MM , J ) , J * l, 13) 

00 435 I = 1, MM 

IF(ENCII) - ENOTIM) 431, 431, 435 
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A3 1 START! 1 1 = START! PM ) 

END ( I t = END ( MM ) 

TYPE I l ) = TYPE IMP) 

NOMST(I) = NCMST(PM) 

STATNMI 1) = STATNP! PM) 

NCEGU) = NCEG(MM) 

STOP l T < I) = STOP I T ( PM ) 

WTOB ( I ) = WTOB(MP) 

DC A32 J = 1« 13 

A32 A ( I , J ) = A ( PP » J ) 

GO TO A 15 

A35 CONTINUE 

A90 I P1NOREAD - 12 * NOCBS) ) 251, 2iO, 210 

251 DO 252 III = l, NCOBS 

252 NST AT ( I 1 1 ) = NOMST1 III) 

210 IP(INDIC) 253, 215, 215 

215 00 250 JJJ = l, NOOBS 
DO 216 III = 1, NOOBS 

IF(NGMST(im.EC.NSTAT< JJJ)) GO TO 220 

216 CONTINUE 

220 1 F ( III. EQ. JJJ) GO TO 25C 
TEMP20 = START I J J J ) 

TEMP21 = ENO(JJJ) 

TEMP22 = TYPE! JJJ) 

TEMP23 = NOMSTIJJJ) 

TEMP2A = STATNP(JJJ) 

TEMP25 = NOEG(JJJ) 

TEMP26 = STOP IT ( JJJ) 

TEMP27 = WT0B1JJJ) 

DO 270 KKK = l, 13 

270 STORIKKK) = AtJJJ, KKK) 

START I JJJ) = START! II I) 

END ! JJJ) = ENC! Ill) 

TYPEIJJJ) = TYPE! II I) 

NOP ST ! J J J ) = NOPS T ! 1 1 I) 

STATNM1JJJ) = STATNP! II 1) 

NCEG! JJJ) = NCEG! Ill) 

STDFIT! JJJ) = STOFIT! II I) 

WTOB ( JJJ ) = WTOB! Ill) 

DO 277 KKK = 1, 13 
277 A ! J J J , KKK) = AIIII, KKK) 

START! Ill) = TEMP2C 
ENC! Ill) = TEMP21 
TYPE! IIII = TEMP22 
NOMST ( III) = TEMP23 
STATNMI III) = TEMP2A 
NCEG IIII) = TEMP25 
STDFIT! Ill) = TEMP26 
WTOB! Ill) = TEMP27 
DC 276 KKK = 1, 13 
276 AIIII, KKK) = STOR! KKK) 

250 CONTINUE 

DO 260 I = 1, NK 

260 WRITEI6.5001) I, STARMI), ENC(I), NOMSTII), TYPE ( I ) , NDEG(I), 
1 1 A! I.J) , J = l,, 13) 

IFIINCCT) 139, 139, 253 
253 T = BEGTIP 

M = 1 

1A N = NOOBS 

IF(NOPRST) 25, 15, 15 
15 IFINCMST(N) - NOPRST) 2C, 25, 20 
20 N = N - 1 

IF1N) 30, 30, 15 

30 WRI T E < 6 , 50 02 ) 
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5002 FORMATdHO, 50X, 33F NO TWO WAY RANGING STATION FOUNO//) 

GC TC 9C00 

25 NUMST = NOMST ( N ) 

NAMBPR(N) = NAMBI NUMST) 

TEMP3 = NAMB(NUMST) * TEMP2 
TO = T ♦ .5 * TEMP2 

NUMST = NOMST(N) 

STNCEL = ST AOEL 1 NUMST ) 

KOUNT = 1 

N2WAY * N 

GO TO 80 

31 CALL CHEBYINC, AC, TO, Y, XST, XEND) 

Y = Y * 1.0D-C6 

Y2 = Y 

YSAVE = Y 

TEMPI = .5 * Y / B640C. 

OELCOOI M,N ) = TEMPI 
TP = TO - .5 * TEMF3 ♦ TEMPI 

I F ( DAB S ( T / TP - l.C) - EPSLGNI 40, 40, 35 
35 TO = TD - TEMPI 

GO TC 31 

40 OELAY(N) = V ♦ NAMB(NUPST) * OELA - STNCEL 

DEL2W = Y ♦ NAMB(NUMST) * DELA - STNDEL 

OELCUT ( M ,N ) = Y - STNCEL 

D2WAY = Y - STNOEL 
WR I T E ( 6 , NAM 1 ) 

T = TP 

I F ( NCPR ST I 42, 41, 41 

42 IF(NOPTS) 26, 26, 65 

26 N = N - 1 

I F( N ) 65, 65, 25 

41 N = NCCBS 

43 I F ( NOMS T ( N ) - NOPRST) 5C, 45, 50 

45 N = N - 1 

I F ( N J 65, 65, 43 

50 TEMP3 = 0.5 * TEMP2 * < NAM8 (NOPRST ) ♦ NAMB ( NUMST ) ) 

TO = T ♦ TEMP3 

NUMST = NCKST(N) 

NAMBPR(N) = NAM8( NUMST) 

STNCEL = ST AOEL < NUMST ) 

KOUNT = 2 

80 NSA V = N 

IF(KCPTS) 77, 77, 76 

76 N =1 

77 NO = NOEG(N) 

XST = START (N ) 

XENC = ENO(N) 

KK = NO ♦ 1 

DO 81 K = 1, KK 

81 AC ( K ) = A (N,K ) 

N = NSAV 

GO TO (31, 51), KOUNT 

51 CALL CHEBY ( NO , AC, TO, Y, XST, XEND) 

Y = Y * 1.0C-C6 

Y4 = Y 

YSAVE = Y 

TEMPI = .5 * Y / G640C. 

OELCOD ( M,N ) = TEMPI 
I F ( Y4 - Y2 ) 52, 52, 55 

52 TO = TO ♦ .5 * TEMF2 

53 CALL CHEBYIND, AC, TO, Y, XST, XEND) 

Y = Y * 1 .00-06 

YSAVE = Y 

Y4 = Y 

DEL AY ( N ) = Y ♦ (NAMBI NUMST) + NAMB( NOPRST ) ) * DELA - STNDEL ♦ OELA 
TEMPI = .5 ♦ Y / 8640C . 

OELOCOl M ,N ) = TEMPI 
WR I TE ( 6 ,NAM3 ) 

I F ( Y4 - Y2 ) 54, 54, 55 
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54 TP = TO - TEMP3 ♦ TEMPI - .5 * TEMP2 

CO TO 56 

55 TP = TO - TEMP3 ♦ TEMPI 

DELAY ( N ) = Y ♦ (NAM8( NUMST ) ♦ NAMB I NCPRST ) ) * OELA - STNCEL 

56 I F ( C A8 S ( T /TP - l.C) - EPSLON) 6C, 60, 57 

57 TO = TC - TEMPI 

GO TC 51 

60 OELAY(N) = OELAY(N) - DEL 2W 

DELCUT ( M ,N ) = Y - 02WAY - STNCEL 
N = N - 1 

IF(N) 65, 65, 43 

65 IF(NCPTS) 82, 82, 83 

82 OC 66 N = 1, NOCBS 

R(N) = OELAY(N) ♦ C / (2.0 * REARTH) 

RTVIIM»N) = R ( N ) * REARTH * 100C.0 

TVT I ME ( M,N ) = T 

66 ROBS(N) = R(N1 * R( N) 

GO TO 85 

83 R ( J J ) = DELAY (N I ♦ C / (2.0 * REARTH) 

RTVKM.N) = R(JJ) * REARTH * 100C.0 
TVT IME ( M,N ) = T 

85 IF(NCPTS) 90, 90, 145 

90 1 F ( NCCB S - 3) 145, 91, 51 

91 I =1 

NLSCIT = 1 

WT< 1) = - 1.0 

CALL RE ADC L ( SYR E , SCNTH , SOAYY, SHOURS, SUETS, SSECCN, KK, T,SE) 

KYR = SYRE 

KMGN = SONTH 

KDAY = SOAYY 

KHR = SHOURS 

KMIN = SUETS 

SECS = SSECON 

Y = YP 


202 

IF(NCMSTII) - 47) 

68, 

67, 

68 

67 

NOSTN = 2 





GO TO 75 




68 

I F ( NCHST ( I ) - 58) 

70, 

69, 

70 

69 

NOSTN = 3 





GO TO 75 




70 

IF(NCMSTd) - 66) 

72, 

71, 

72 

71 

NOSTN = 4 





GO TO 75 




72 

I F ( NCMS T ( I ) - 68) 

74, 

73, 

74 

73 

NOSTN = 1 





GO TO 75 

74 MR 1 T E ( 6 , 5003) NOMSTII) 

5C03 FORMAT ( 1H0 , 11HSTATI0N N0..2X, 13, 2X, 58HIS NOT INCLUDED IN THE 
1STAT ION CONSTANTS, PROGRAM HALT EC. // 1 8X , 65HCHECK THE STATION NO., 
20R ACC NEW STATION TO PROGRAM ANC TRY AGAIN) 

GO TC 9C00. 

75 MDAY = KDAY 

I F ( KMON - IMON) ICO, 11C, ICO 
100 MDAY = MC AY ♦ NDAY(IMON) 

110 LROh = MDAY - IDAY ♦ 1 
ALAMCO = ALAMDM(LRCW) 

ASECS = 3600 * KHR ♦ 60 * KM I N 
ASECS = ASECS ♦ SECS 
CELT = ASECS / 360C.C+C0 
ALAMCE = S T ACON ( NO STN, 1) 

THETAC = S T ACON ( NOSTN , 2 ) 

ALT = STACOMNOSTN, 3) 

SNTHEC = OS IN ( THE T AD ) 

CSTHEC * DCCS(THETAD) 

CAPC = CSTHED ** 2 ♦ 1(1. - F) * SNTHED ) ** 2 
CAPC = l.C / DSQRT ( C APC ) 
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CAPS = CAPC * (1.0 - F) ** 2 

ALT = ALT * 4.7 7865C-C8 

TEMP = (CAPS ♦ ALT) / (CAPC ♦ ALT) 

TEMP = TEMP * SNTHED / CSTHED 
THETAG = OATAN (TEMP) 

SNTHEG = OS IN ( THETAG ) 

CSTHEG * DCOS( THETAG) 

DELTA =ALAMDO ♦ ALAMOE ♦ OELT * 2. 625 161 333D-0 1 
SNDELT = OS IN( OELT A) 

CSOELT = OCOS(OELTA) 

RHOCAP = (CAPS ♦ ALT) ♦* 2 

RHOCAP = RHOCAP ♦ SNTHEO ** 2 ♦ {(CAPC ♦ ALT) * CSTHED) ** 2 
RHOCAP = DSQRT (RHOCAP) 

COMPUTATION OF INERTIAL CCORDINATES OF OBSERVATION POINT 

XM = RHOCAP ♦ CSTHEG ♦ CSOELT 

VM = RHOCAP ♦ CSTHEG * SNDELT 

ZM = RHOCAP * SNTHEG 

XT( I) = XM 

YT< I ) = YM 

ZT( I) = ZM 

I =1*1 

I F ( I - NOOBS) 202, 202, 116 
116 I =1 

I F ( NO I T ER - 1) 118, 118, 400 
118 CALC COMP < XT, YT, ZT, RCBS, XOF, YOF , ZOF) 

X = XCF 

Y = YOF 

Z = ZOF 

MR I TE ( 6 , 5000 F, TOL, DELA, NOOBS, NOPRST, OELMIN, REARTH, C 
5000 FORMAT ( 1H1 , 50X, 25H ATS EXPERIMENTAL RANG I NG//54X , 

1 16H INPUT PARAMETERS //4X, 20HE ARTF • S FLATTENING = E15.8, 5X, 
223HC0NVERGENCE TOLERANCE = E15.8, 5X, 6H0ELA = E15.8, 

32X, 3HSEC//4X, 27HN0. OF OBSERVING STATIONS = 15, 5X, 

422HNC . CF PRIME STATION = 15, 5X, 8HDELTAT = E15.8, 2X, 4HMINS// 
518H RAO I US OF EARTH = E20.8, 2X, 3HKM. , 5X , 16HSPEEC OF LIGHT = 
6E20.8, 2X, 10HKM. / SEC.//) 

WRITE(6,5010) X, Y, Z 
5C10 FORMAT ( 1H0, 

166X, 3HX = E15.8//1CX, 

259H IN IT I AL ESTIMATE OF INERTIAL COORDINATES OF SATELLITE Y = 

3 E15.8//66X , 3HZ = E15.8) 

MR I TE ( 6 , 5005) 

5005 FORMAT ( 1H 1 , 50X, 6HCUTPLT // 

1 9X, 4HTIME, 26X, 1 9HSU E-SATELL I T E POINT//IX, 

222HYR MO DAY HR MIN SEC, 6X, 13HLCNG I TUOE ( +E ) , 2X, 
324HLATITUDE(+N) HE IGHT (KM) // ) 

400 CLMCKI, 1) = X - XT( I) 

CLMCKI ,2) = Y - YT( I) 

CLMOK 1,3) = Z - ZT( I) 

RCOMP = CLMQI (1,1) ** 2 * CCMQI(I,2) *♦ 2 ♦ CLMCKI, 3) ♦* 2 

RC2 = DSCRT ( RCOMP I 

OELR ( 1 , 1) = .5 * ( ROBS ( I ) - RCOMP) 

I F ( I - NOOBS) 120, 130, 130 
120 1=1+1 
GO TO 4C0 

130 CALL GLSP(CLMQI , NOCBS, 3, OELR, 1, DLTAQI, RESIO, SUM, MT, 0, 
1ST0ERR) 


X 

* 

X 

♦ 

OLTAQK 1,1) 

Y 


Y 

■f 

DLTAQK2, 1) 

Z 

= . 

Z 

♦ 

OLTAQK 3, 1) 

xs 

= 

X 

* 

X 

YS 

= 

Y 

* 

Y 

ZS 

= 

Z 


Z 

NLSQIT 

= 

NLSQIT ♦ 1 

SUM 

= 

0 

• 


00 3C0 

I 

= 

1 

, NOCBS 
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300 SUM = SUM ♦ <2.0 * CELR(I,1I) ** 2 

FIT = OSQRT ( SUM / NCOBS) 

I = l 

I F ( NL SQ l T - 20) 350, 35C, 8C00 
350 IFIFIT - TOL) 500, 500, ACO 

BEGIN COMPUTATION OF SUE-SATELLITE POINT 

500 RHOS = XS ♦ VS + ZS 

I =0.0 

OTHEG = 0.000 

TES = F * (2.00C - F) 

T2 = ( 1. - F) ♦* 2 

RHO = OSQRT IRHOS) 

SNTFEG = Z / RHO 

CSTHEG = OSCRT ( ( XS + YS ) / RHOS) 

TEMPI = RHC * CSTHEG 

CSDELT = X / TEMPI 

SNDELT = Y / TEMPI 

OELTA = 0ATAN2I SNDELT, CSOELT) 

XLONG = OELTA - ALAKDO -.26251613330+00 * CELT 
DEL = C AT AN ( Z/OSQRT (XS ♦ YSI) 

501 THETAG = DEL - OTHEG 
1=1 + 1 

T 10 = OCOS(THETAG) 

Til = T10 ♦ T 10 

RC = OSQRT ( ( 1. COC - TES) / (1.000 - TES * Til)) 

XL AT = DATAN(DTAN(THETAG) / T2) 

T12 = XL AT - THETAG 

T 1 3 = DSINIT12) 

T 14 = T 1 3 ♦ T 13 

T 15 = OCOS ( T 1 2 1 

HT = OSCRT (RHO S - RC * RC * T1A) - RC ♦ T15 

T 20 = HT * T 1 3/ RHO 

OTHEG = D ARS I N ( T2 0 ) 

IF( I.LT.10) GO TO 5CI 
XLAT = XL AT * RACOEG 

XLONG = XLONG * RADCEC 

HT = HT * AE 

520 I F ( OABS ( XLONG) - 18C.0) 510, 510, 700 

700 IF ( XLONG ) 701, 702. 702 

701 XLONG = XLCNG ♦ 360.0 
GO TO 510 

702 XLONG = XLONG -36C.0 

510 WR I T E ( 6 ,5004 ) KYR, KMCN , KOAY , KFR , KM IN , SECS, XLONG, XLAT, HT 
5C0A FORMAT! 1H0, 12, IX, 12, 2X , 12, IX, 12, 2X, 12, 2X, F5.2, 5X, F8.3, 
17X , F7 o 3 , AX, F 10.2 ) 

WRITE(6»9035) X, Y, Z 

9C35 FORMAT! 1H0, AH X = E20.E//AH Y = E20.8//AH Z = E20.8) 

1A5 T = T + DELTA! 

NOITER = NO I TER + 1 

YP = Y 

M = M ♦ 1 

JJ = JJ ♦ 1 

IF(T - ENOTIM) 150, 150, 9000 
150 IF(NOPTS) 1A, 1A, 25 
9000 IF(NCPTS) 170, 170, 160 
160 N = N ♦ 1 

JJ = 1 

IF(N - NOOBS) 5, 5, 165 
165 N =1 

170 NTOT = M - 1 

WRITE(9,1000) NTOT 
WRITE (6,1000)NTOT 

1C00 FORMAT! 15) 9 

M =0 
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8999 H * M ♦ 1 

WRIT£(9,1001XRTV1IM,NN), NN => 1, NQCJBS) 
WR1TE(6,1001)(RTVMM,NN), NN = X, NOCBS) 
WRITE(9,1001)(TVTIME(M,NN), NN * I, NOOBS) 
WRITE(6,1001)(TVTIM£(M,NN), NN = 1, NCOBS) 
1C01 FORMAT! 3025.15} 

IF(M - NTOT ) 8999, 9010, 9010 
9C10 00 9C50 M » 1, NTCT 
DO 9C50 L = 1, NOCBS 

CONVERT RANGE MEASUREMENTS TO OUL UNITS 


IOBN = 10 BN ♦ l 

ROUT = RTVl(M.L) * 1 .0C-07 

TIME = TVTIME(M.L) 

T 30 = (NAMBPR(L) * DELA) / (2.0DO * 86900. ODO) 

TIME = TIME ♦ T3C ♦ CELDOOIM.L) 

IF< L.NE.N2WAY) TIME = TIME ♦ T30 - RT V l ( M, N2 W A Y ) / (1000.0 ♦ C * 

1 869C0.0) 

CALL READCUSYRE, SCNTH , SO AY Y , SHOURS, SUETS, SSECON, KK, TIME, 
1SE J 

KYK = SYRE 

K MON = SCNTH 

KOAY = SO AY Y 

KHR = SHOURS 

KM I N = SUETS 

SECS = SSECON 

IF(CFL0AT(KYR/9) - CFLOAT(KYR) / 9.000) 99, 92, 99 
92 NCAYI2) = 29 

99 NOOAY = 0 

LL = KMON - 1 

00 95 11= 1, LL 
95 NCOAY = NCCAY ♦ NCAY(II) 

NOOAY = NODAY ♦ KCAY 

NSJCS = (KYR - 58) ♦ 365 ♦ (<KYR - 1) / 9 - 19) ♦ NCOAY ♦ 109 

SJOS = ICO. 00 * (OFLCAT ( NS JOS ) ♦ UFLOAT ( KHR ) / 29.00 
1 ♦ CFLCAT(KMIN) / 1990.00 ♦ SECS / 86900.00) 

OBC V = 0 • CO 

RCBVT = HTCB(L) 

STNAME = STATNM(L) 

SOFT = STOFIT(L) 

9050 WRITE (29ISJCS, STNAME, RCUT , TR ANCM, F I ORP ,OBC V , I SA T , I OBN , ROBWT ,COC , 
1P0C ,SOFT,FNRRES,TMER,TCPV,IESN, ICRPN, IOT ,IEF, I T T I , I TC I , IOC I , ICC I 
GO TO 5 

8C00 WRITE(6,3050) 

3C50 FCRM AT ( IHO , 95X , 31E PRCGRAM HALTED, NO CONVERGENCE) 

8S98 REWIND 29 

WRITE (6,1300) 

1300 FORMAT ( • 1 OU T • , 19X , • S TN AME OBSV* , 16X, 'FEEDBACK CORR VAL', 

1 7X , ' SA T NO, 08 NO, OB WT',8X,'ST. DEV. CBTP ECFG CC I • ) 

CO 1302 J=l, 32000 

RE AC (29,ENC=5210)S JOS, STNAME, ROUT, TRANOM.F I ORP, OBC V, I SAT, l OBN, 
10BWT,C0C,P0C,SDFT,FNRRES,TMER,TCRV,IESN, IORPN, 1CT, I E F , 1 T T I , I T C l , 

2 IOC I 

1302 WRITE (6, 1301 IS JOS, STNAME, ROUT, I 1 1 ( 1 ) , 1 1 1 ( 2 ) , OBC V , I SAT, IOBN, 
10BWT,SDFT,ICT,IEF,ICCI . 

1301 FORMAT (IX, 020. 12, IX, A8, IX, 020. 12, IX, 215, IX, D19. 6, IX, 17, 16, 

11X,C12.9,1X,C12.9,IX, 12 ,1X, 12, IX, 16) 

5210 STOP 
END 

/♦ 


0583 CAROS 
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SUBROUTINE COMP 

SUBROUTINE COMPIXT, YT, ZT, ROBS, XOF, YOF ■ ZOF) 

IMPLICIT REAL*8!A--I,3-Z,$) 

REAL*8 LAMBDA 

DIMENSION XI 3) , YT 31 ,2(3 ) , XO I 2 ) , YOI 2 ) ,Z0(2) ,RO(3,2) ,R02(2I .CHECK! 2) 

1, XT(IO), YT 110), ZTI10J, ROBSIIOI, R(3,3), RH0I3), R2I3), 

2 RHON (3,3) 

DO 25 1=1,3 
Ml, I) = -XT! I I 
R ( 2 , I ) = - YT! I ) 

R( 3, I ) = - ZT ( I ) 

RHO! I I = DSQRT (ROBS! I ) ) 

XI I )=R( 1, I ) 

Y( I )=R< 2, I ) 

25 Z( I )=R( 3, I I 
DO 35 1=1,3 

35 R2( I ) = DOT ( R ( 1, I ),R( 1, I ) ) 

E2l=.5D0«( RHO! 2 )**2-RH0(l )**2-( R2(2)-R2t l>)) 
E31=.5D0*(RH0(3)**2-RH0(1)**2-1R2(3>-R2(1)>» 

DELTAl= ( Z< 3I-Z1 1)1* (Y(2 )-Y( II )-(Z(2)-Z( 1)) *( YI31-YI 1)) 
A=((X(2)-X(in*(Y(3)-Y(m-!XI3)-X (!))*! Y( 2) -Y ( 1 * ) ) /0ELTA1 
B*I E31*( Y! 2 ) -Y C 1) ) - E2 1* I Y( 3 ) -Y( 1 I ) J/DELTAl 
DELTA2=-DELTA1 

C = (!X(2)-X(in*!Z(3)-Z(tn-(X(3)-X(tn*!ZI2)-Z(l)) I/DELTA2 
D= ( E 31 * ( Z ( 2J-ZI 1) )-E21*!Z(3)-Z< 1 ) ) 1/DELTA2 
EPSl=A*A+C*C+1.00 

EPS2 =2.D0*IA*B»C*0«lllK»rtl!*Wllll 
EPS3 =B*B»D*D+2.30*D*Y( 1 I +2. DO*B*Z (1) +R2 111 -RHO(l) **2 
RAD =EPS2 **2-<,.D0*EPSl*EPS3 
XO! 1 ) = (-EPS2+DSQRT( RADI )/<2.00*EPSl> 

XO! 2) = ( -EPS2-DSQRT(RAD) 1/(2.D0»EPS1 I 
DO 30 1=1,2 
Y0< I >=C*XO( I I +D 
Z0( I)=A*XO( I ) +B 
RO ( 1,I)=X0( II 
R0( 2, I J =Y0( I ) 

RO(3»I)=ZO( I) 

XOF = XO( II 
YOF = YO ( 1 1 
ZOF = ZO( I ) 

DO 75 X = 1 » 3 
DO 75 J = 1 ,'3 

75 RH0N!J,K1=R( J,KHRO!J,l ) 

I F ( DOT ( RH0N,R) /( DS3RT (R2) *RHO( 1) ) .LE.O.DOt RETURN 
30 CONTINUE 
RETURN 
ENO 

CDOTT 

FUNCTION DOT ( A , B ) 

IMPLICIT REAL*8(A-H,0-Z,$) 

DIMENSION A ( 3) , B ( 3) 

D0T=0. DO 
DO 1 1=1,3 
1 DOT = DOT + A( I )*B( I ) 

RETURN 

ENO 0058 CARDS 
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c 

c 

c 


SUBROUTINE CHEBY 


■‘• 00 

TO= I .00+0 

Uc(n*TO»ciz>*Ti 

IF t NO* LT. 2 ) RETURN 

n t=no+i 
00 10 

TN=2.0*X*T1-T0 

V = Y*C(<1* tn| 

to = ti 
13 T1*TN 
RETURN 
P MO 


3319 CURDS 


SUBROUTINE DOT 


l 


IMPLICIT r 

0 1 MENS I ON 
OOT-0.00 


E&L*9 1 A-R 
M 3) .BO) 


00 V l= l 
00T = 00 T «■ 




3 

( l )* 3 ( I) 


return 

ENO 

//OBJECT DO * 
/* 

//OBJECT DO * 


, 0 -Z .*> 


0013 C&RDS 


C 

C 

C 


M ‘“ l t , s< «,s OT ™. S o.vv.s-«s.s U e.s.ss E co N . 

SUBROUTINE REAQCL 

® g0&me mses8 ak 

<,*.6065. 00/ 

S UNA=SSUNA 
E = + l .DO 
1 = 1 

„ s,ffiiusriSi‘”- K J,28 ’„ J0 


66 


FCPT0010 


FCPT0020 

FCPT0030 

fcptooao 

FCPT00T0 
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67 IF (YJ(28)-S UN A)12« 68*68 

68 1=2 

71 IF < YJ< I )-SUNA+.5D0 169, 70,70 

69 1=1+1 

GO TO 71 

70 YRE=YM1I-1) 

100 S>UNM=SUNA-YJ( 1-1) 

KSUNM=SUNM 
SUNM 1 = K SUNM 
IYRE*YRE+.01D0 
I F ( MODI I YRE , 4) ) 77, 78,77 
78 DM ( 2 ) =2 9. DO 

GO TO 105 
77 OM ( 2 ) = 2 8« DO 

1 05 1=1 

SM=OM( I 1 + 1. DO 
73 IF (SM-SUNM)8l,81,72 
81 1=1+1 

IF (1-13)210,72,213 
210 SM= SM + OM( I ) 

GO TO 73 

72 SD=SUNM-( SM-OM( I) J+l.DO 
ONTH= I 

201 KSD= SD 

DAC = SD 

OAYY=KSO 

SDAC=DAC 

SDAC= I D I NT ( SDAC ) 

DO AC=SD AC 
SD=OAC- DDAC 
IF (SD 1220, 221. 221 

220 SD= 1 .DO+SD 

221 H0UR=$D*2A.00 
KHOUR=HOUR 
HOURS=KHOUR 

UE= (HOUR-HOURS 1*60.00 

KUE =UE 

UETS=KUE 

SECON=(UE-UETS)*60.DO 
IF (ONTH)207,208,207 
208 ONTH= 1 2 . DO 
207 CONTINUE 
99 SYRE=YRE 
SONTH=ONTH 
SDAYY=D AYY 
SHOURS=HOURS 
SUETS=UETS 
SSECON= SECON 
SSUNA=SUNA 
SE = E 
RETURN 
END 


FCP70080 
FCP70090 
FCP70100 
FCP701 10 

FCP73130 

FCP701A0 

FCP70150 

FCP73160 

FCP70170 

FCP70190 

FCP70210 

FCP70220 

FCP70230 

FCP702A0 

FCP70250 

FCP70260 

FCP70270 

FCP70290 

FCP70300 

FCP70310 

FCP70320 


FCP703A0 

FCP70350 

FCPT0360 

FCP70370 

FCP70380 

FCP70390 

FCP70A00 

FCP73A10 

FCP70A20 

FCP70A30 

FCP70AA0 


FCP71A30 

0080 CARDS 


SUBROUTINE GLSP 


SUBROUTINE GLSP ( A , MM, NN , 3 , I PP , X , U, SUM, WT , I NVRS , STDERR ) 005 

IMPLICIT RE AL*8 (A-H.O-Z) 

DIMENSION A ( 10,3), B ( 10,1), X (3, A), U ( 10,1), SUM (A), 

1 W (3, A), WT ( 10), STDERR (3) 

M=MM 008 

N=NN 009 

IP=IPP 010 

I NVRSE= I NVR S Oil 

AX=B, NORMAL EQUATIONS, A IS RECTANGULAR MATRIX, WITH M ROWS.N COLUMNS 012 
WX = W1 „NEW SET OF NORMAL EQUATIONS, W=TR(A)*A ,Wl= TR ( A) *B. . 013 

TO FORM MATRIX W .WITH N ROWS AND N COLUMNS. ( TR= TRANSPOSE ) OlA 
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C WT = -1. MEANS NO WEIHGTS IN PROG. ASSIGN WEIGHTS = 1. 015 

C STOERR = STANDARD ERRORS OF UNKNOWNS ( X ) . 016 

C INVRSE = 1. MEANS STO. ERRORS ARE CALCULATED. 017 

C INVRSE = 0. MEANS NO CALCULATIONS OF STD. ERRORS OF ( X )♦#**#**** 018 

KK = N + 1 019 

IFIM-N) 80,83,81 020 

81 DO 5 J = l , N 021 

DO 5 J J = l , N 022 

WIJ.JJ) = 0. 023 

IF ( WT ( 1 ) ) 611, 622, 622 

C NO WEIGHTS CASE 025 

611 00 51 1= l,M 026 

51 WIJ.JJ) = WIJ.JJ) ♦ A I I , J ) * AII.JJ) 027 

GO TO 5 028 

C WEIGHTS ARE PRESENT 029 

622 DO 52 1=1, M 030 

52 WIJ.JJ) = WIJ.JJ) ♦ All . J ) * AII.JJ) * WTIU 031 

5 CONTINUE 032 

C TO FORM MATRIX W1 , WHICH IS STORED IN N* 1 TO N*P COLUMNS OF W. 033 

C W1 .WITH N ROWS AND IP COLUMNS. 035 

IF I INVRSE ) 201,203,201 0035 

203 KIP = IP 036 

GO TO 200 037 

201 KIP = 1 038 

200 DO 8 K = 1, KIP 039 

DO 7 J = 1 , N 050 

WIJ.KK) = 0. 051 

IF I WT I D) 612, 623, 623 

612 DO 511 1=1, M 053 

511 WIJ.KK) = W<J,KK)+AII,J J* BII.K) 055 

GO TO 7 055 

623 DO 522 1= 1,M 056 

522 WIJ.KK) = W< J,KK)+A(I,J )* B(I,K) * WTI I ) 057 

7 CONTINUE 058 

8 KK = KK ♦ l 059 

GO TO 106 050 

83 DO 12 1=1, N 051 

IF (WTI 1)1 613, 625, 625 

613 DO 512 J=l,N 053 

512 W( I , J) = A( I , J) 055 

GO TO 12 055 

625 DO 523 J=1,N 056 

523 W( I , J) = A< I ,J) * WTI I ) 057 

12 CONTINUE 058 

KKK = N+IP 059 

J=1 060 

IF I INVRSE ) 501,503,501 061 

503 KKK = KKK 062 

GO TO 500 063 

501 KKK = KK 065 

500 DO 13 JJ= KK.KKK 065 

IF (WTI 1) ) 615, 625, 625 

615 DO 513 1= 1 , N 0067 

513 WI I ,JJ) = B I I » J ) 068 

GO TO 13 069 

625 DO 525 1= l,N 070 

525 WII.JJ) = BII.J) * WT ( I ) 071 

13 J=J+1 072 

106 IF I INVRSE ) 71,72,71 0.73 

71 K2 = N *• 2 075 

KIP = N ♦ IP 075 

KK = N ♦ 1 076 

DO 75 I =1 , N 077 

DO 75 J = K2.KIP 078 

IF I J - I - KK ) 73,75,73 079 


39 



1 


75 

73 

75 

72 

C 

c 

c 


82 

80 

103 

C 

c 

c 

105 


15 

16 


18 

722 

723 


725 

105 


HU.JI = 1. 

30 TO 75 
H(I . JJ = 0. 

CONTINUE 

CALL TRIANG (W.N.IP, 
WHERE W IS THE MATRIX 


DET 1 

WITH N ROWS AND ( N+ I P 1 COLUMNS 


R I ZED 

TO TEST THE SINGULARITY OF MATRIX W.. W = TR t A1 *A 
IF ( DET I 82,80,82 


TO BE TRIANGULA 
, W*X = wl. 


080 

081 

082 

083 

085 

085 

086 

087 

088 


CONTINUE 
30 TO 103 
SUM ( 1 ) = -l. 

GO TO 105 

CALL SOLVE (W,N,IP,X) 

WHERE W IS THE TR I AN3UL AR l ZEO , 

N ROWS AND IP COLUMNS. 

TO COMPUTE RESIDUAL MATRIX , U « A*X 
DO 16 1=1, M 
DO 16 K = 1 , I P 
UU.KI =0. 

DO 15 J = 1 , N 

UU.KI = UU.KI ♦ AII.JJ * X ( J , K ) 
UU.KI = B( I ,K J - UU.KI 
DO 18 K = l, IP 
SUM ( K I = 0. 

DO 18 1=1, M 

SUM ( K I =SUM!K) ♦ ( U( I, K) 1**2 
IF ( INVRSE 1 722,105,722 
00 723 I = l.N 
DO 723 J» 2, IP 
XU.Jl = XU.Jl * 

00 725 I = 1 , N 
STOERR ( 11=0. 

II = I ♦ l 
STDERR I 1 1 = OSORTI 
RETURN 
END 


089 

090 

092 

093 

MATRIX WITH 095 

095 

- B , U IS RESIDUAL MATRX t M, P ) 096 

097 

098 

099 

100 
101 
102 
103 
105 

105 

106 

107 

108 
109 

111 
112 
113 


0116 CARDS 


X IS THE SOLUTION 


( SUM (II / DFLOAT ( M-Nl 
XU, II 11 


C 

C SUBROUTINE TRIANG 


SUBROUTINE TR I ANG ( A ,M, N ,DET 1 120 

IMPLICIT RE AL*8 (A-H.O-ZJ 
DIMENSION A (3,51 

CTRIAN3 TO BE USED WITH SUBROUTINE GLSP H q 

C . WHERE A IS THE MATRIX WITH M ROWS AND M*N COLUMNS TO 121 

C BE TRIANGULARIZED. 122 

DET=l. tz * 

MMl=M-l 125 

MPN=M«-N 126 

DO 5 J = l , M M 1 127 

MAXX = J 128 

VALUE = DABS ( A(J,J!1 

JPl=J+l 1 30 

DO 1 K= JP1 » M 151 

IF (VALUE - DABS ( A( K, J) 11 2, l, l 
2 VALUE = OABSI A(K,J11 

M AX X=K 155 

1 CONTINUE 155 

DO 3 L = J, MPN 136 

TEMP = At MAXX, L 1 137 

A ( MAXX , L 1 = A ( J , L 1 138 

3 A ( J , L 1 = TEMP 139 

IF t MAXX- J 1 7,6,7 1*0 

7 DET = - DET 1*1 


40 



ooooo noo 


6 OET = OET*A(J,J) 1« 

ROW = A ( J , J 1 143 

00 4 L'J.MPN 144 

4 A(J,U) = A<J,L)/R0W 145 

DO 5 K=JPl,M 146 

ROW = AIK.JI 147 

00 5 L s J » HPM 148 

5 A ( K , L I = A(K,L1-R0W*A<J,L) 149 

OET = DET*A(M,M) 150 

ROW * AIM, Ml 151 

DO 9 L = M , MPN 152 

9 A I M , L ) = A(M,L)/R0W 

RETURN 154 

EN0 0041 CARDS 


SUBROUTINE SOLVE 


SUBROUTINE SOL VE I A, M, N, X) 159 

IMPLICIT RE AL*8 (A-H.O-Z) 

DIMENSION A (3,4), X I 3,41 

TO BE JSEO WITH SUBROUTINE GLSP 158 

WH£RE A IS THE TRIANGULAR! ZED MATRIX WITH M ROWS AND 160 

M*N COLUMNS. 161 

X IS THE SOLUTION MATRIX WITH M ROWS AND N 162 

COLUMNS. 163 

MM1=M-1 165 

DO 2 L=l.N 166 

MPL=M+L 167 

DO 2 K = l , MM 1 168 

MMK=M-< 169 

MMKPl=MMK*l 170 

X(M,LI=A(M,MPL)/A(M,M) 171 

SUM=0. 172 

DO 3 I = MMK P 1 , M 173 

3 SUM=SUM*A( MMK, I )*X( I,L) 174 

2 X I MMK, L ) =A( MMK, MPL) -SUM 175 

RETURN 176 

END 177 


0024 CARDS 


41 



